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Abstract. In this paper, discrete analogues of Euler-Poincare and Lie-Poisson 
reduction theory are developed for systems on finite dimensional Lie groups G 
with Lagrangians L : TG —> K that are G-invariant. These discrete equations 
provide "reduced" numerical algorithms which manifestly preserve the sym- 
plectic structure. The manifold G X G is used as an approximation of TG, and 
a discrete Langragian L : G X G — > M is construced in such a way that the G- 
invariance property is preserved. Reduction by G results in new "variational" 
principle for the reduced Lagrangian I : G — * M, and provides the discrete 
Euler-Poincare (DEP) equations. Reconstruction of these eq u ations r ecovers 
the discrete Euler-Lagrange equations developed in [ MPS 9a , WM 97] which 
are naturally symplectic-momentum algorithms. Furthermore, the solution of 
the DEP algorithm immediately leads to a discrete Lie-Poisson (DLP) algo- 
rithm. It is shown that when G = SO(n), the DEP and DLP algorithms for 
a particular choice of the discrete Lagrangian L are equivalent to the Moser- 
Veselov scheme for the generalized rigid body. 
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1. Introduction 

The goal of this paper is to develop structure preserving numerical integrators 
on the reduced space of a mechanical system whose configuration space is a Lie 
group G, and whose Lagrangian L : TG — > K. is either left or right invariant by the 
group action. In particular, we shall develop the discrete analogue of Euler-Poincare 
theory by following the variational approach introduced by Marsden, Patrick, and 
Shkoller |MPS 98] for the construction of discrete Euler-Lagrange equations that 
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naturally preserve the symplectic structure and the momentum mappings of the 
Lagrangian system. 



In our setting, the results of [MPS 95] may be described as follows. Given a 



Lagrangian L : TG —> R, form the action S on curves g : [a, b] — > G defined in a 
chart by 

S(g(t))= f L(g*(t),g*(t))dt. 

J a 

Allowing for arbitrary variations 8g, not constrained to vanish on {a, 6}, a compu- 
tation of the first variation of S leads to 

•«')- />(! -iw)* + P 3 i (L1) 



The last term of (1.1) is a linear pairing of dL/dg l , a function of g 1 and g l , with 
the tangent vector Sg l . Thus, one may consider it to be a 1-form 9l = (dL/dq l )dq l 
on TG, and the symplectic structure is then defined by 

u L = —d6 L . 

Applying the operator d 2 =0 to S, restricted to the space of solutions of the Euler- 
Lagrange equations, shows that the flow F t of the Euler-Lagrange equations con- 
serves the symplectic form; namely, F^cul = ^>l- Next, let g denote the Lie algebra 
of G and define the momentum mapping : TG — > R for each £ 6 g corresponding 
to the tangent lift of the right (or left) action of G on itself by = £tg-I 9l, where 
£tg is the infinitesimal generator of £ 6 q on TG. Then, the variational principle 
together with the infinitesimal invariancc of the action restricted to the space of 



solutions, immediately leads to the fact that F^J^ — J^. See [MPS 98] for details. 

Hence, this variational approach can be used to obtain a symplectic-momentum 
integrator by discretizing TG and forming a discrete action sum. For every choice 
of discretization, a unique discrete symplectic structure is obtained, and the algo- 
rithm given by the discrete Euler-Lagrange equations is guaranteed to preserve this 
structure as well as the momentum mappings associated with it. Our goal is to 
apply the reduction procedure in this discrete setting, restrict the Lagrangian to 
the reduced space, and derive the algorithm which preserves the induced structure. 

Our procedure results in the discrete Euler-Poincare equation which defines an 
algorithm on the reduced space that is shown to be equivalent to the discrete 
Euler-Lagrange equations in the sense of reconstruction. This reduced algorithm 
is used together with the coadjoint action to advance points in g* = T*G/G and 
thus approximate the Lie-Poisson dynamics. In subsequent papers, we shall make 
the extension to the more general setting of Lagrangian reduction of a G-invariant 



system on TQ (see, for example, Cendra, Marsden, and Ratiu [CMR95]), for a 
general manifold Q, as well as to the case of dynamical systems defined on Lie 
algebras. 

2. The discrete Euler-Poincare algorithm 

In this section we develop the discrete Euler-Poincare reduction of a Lagrangian 
system on TG. We approximate TG by G x G and form a discrete Lagrangian 
L:GxG^l from the original Lagrangian L : TG — » R as 

H9k,9k+l) = L(K(g k ,g k+1 ),x(gk,gk+i)), 
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where k and x are functions of {g k ,gk+i) which approximate the current config- 
uration g(t) S G and the corresponding velocity g(t) S T g G, respectively. We 
choose particular discretization schemes so that the discrete Lagrangian L inherits 
the symmetries of the original Lagrangian L: L is G-invariant on G x G whenever 
L is G-invariant on TG. In particular, the induced right (left) lifted action of G 
onto TG corresponds to the diagonal right (left) action of G on G x G. 
Having specified the discrete Lagrangian, we form the action sum 

N-l 

§ = L (Sfc,5fc+i) 

k=0 

and obtain the discrete Euler-Lagrange (DEL) equations 

D 2 h(g k -i,g k ) + D l L(g k ,g k+1 ) = 0, (2.1) 
as well as the discrete symplectic form wl given in coordinates on G x G by 

^ = * in i d & A d 9l+i ( 2 - 2 ) 
a 9k a 9k+i 



by extremizing § : G N+1 — > K with arbitrary variations. It is shown in [ MPS 98| 



that the flow Ft of the DEL equations preserves this discrete symplectic structure. 
We remark here that the original canonical symplectic form oj is also preserved by 
this flow. Indeed, as the discrete Legendre transformations define a local symplec- 
tomorphism, we obtain that u(t) — iTL -1 (u^(t)) — FL,^ 1 ^^^)) = w(0). 

The discrete reduction of a right-invariant system proceeds as follows. The 
induced group action on G x G is simply right multiplication in each component: 

9 ■ (9k,9k+i) >-> (g k g,g k +ig), 
for all g, g^, g k +i G G. Then the quotient map is given by 

ir:GxG^(GxG)/G^G, (g k ,g k+ i) h- g k g£ v (2.3) 
We note that one may alternatively use gk+ig^ 1 instead of gkgk+i'i °ur choice is 



consistent with other literature (see, for example, MPS 98 |). The projection map 
( j2.3| ) defines the reduced discrete Lagrangian i : G — » K for any G-invariant L by 
I o 7r = L, so that 

K9k9kli) = H9k,g k +x), 
and the reduced action sum is given by 

S = Y Kfkk+l), 
k=0 

where fkk+i = gkgk+i denote points in the quotient space. A reduction of the DEL 
equations results in the discrete Euler-Poincare (DEP) equations. We state this 
as the following theorem. 

Theorem 2.1. Let L be a right invariant Lagrangian on G x G, and let £ : (G x 

G)/G = G — > K be the restriction of L to G given by iigig^ 1 ) = ^(91,92)- For 
any integer N > 3, let {(gk, g k +i)}k=o be a sequence in G x G and define fkk+i = 
3fc3fc+i t° be the corresponding sequence in G. Then, the following are equivalent. 
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(1) The sequence {(gk, gk+i)}k=a * s an extremum of the action sum § : G N+1 — > 
R for arbitrary variations Sgk — (d/de)\og k where for each k, e i— > g\ is a 
smooth curve in G such that g^ — gk- 

(2) The se quen ce {{gk, g k +i)}k=o satisfies the discrete Euler- Lagrange equa- 
tions (Jlj). 

(3) The sequence {fkk+i}^=o * s an extremum of the reduced action sum s : G — > 
R with respect to variations dfkk+i, induced by the variations 5gk, and given 
by 

Sfkk+i = TRj^iSgkgk 1 - Ad /fcfc+1 -8g k+1 g^l x ). 

(4) The sequence {fkk+i}k=o satisfies the discrete Euler- Poincare equations 



-£'(f k - lk )Ad fk _ lk TR fk _ lk +e'{f kk+1 )TR fkk+1 = 



(2.4) 



for k = 1, N — 1, where the operators act on variations of the form $ k = 

Sgkgk 1 - 

Proo f. We be gin with the proof that (1) and (2) are equivalent following [ MPS 98| 
and [ WM 97 ] . One computes the first variation of the discrete action § with varia- 
tions that vanish on the set k = {0, N}. Thus, 



d 
~dc 



c=0 



d_ 

de 



N-l 



e=0 fc= 



$>(<M+i) 

k=0 

^2 D iH9k,gk+i)5gk + ^2 D 2 h(gk,gk+i)Sg k +i 

k=0 k=0 
N-l N-l 

D iM9k, gk+i)Sg k + ^2 D 2Hgr-i,g r )Sg r 



N-l 



N-l 



k=l 
N-l 



^2 ( D iH9k, gk+i) + D 2 L(g k -i, gk)) Sg k , 



k=l 

where we have used the discrete analogue of integration by parts which simply 
shifts the sequence gk <— * g r where r = k + 1. Since for each k = 1, ...N — 1, the 
variations 8gk are arbitrary, this establishes the DEL algorithm. We remark that 
choosing variations which do not vanish at k = and k — N defines two 1-forms 
whose exterior derivative is the unique symplectic 2-form given in ( ]2.2| ). 
To see that (1) is equivalent to (3), notice that since L = £ o ir, 



d_ 

de 

Now for (3) (4), we compute 



s(M k+ i)- j e 



de 



N-l 
k=0 



Jk9k+ 



r 1 ) 



and find that 

s{flk+i) = e '(fkk+i) [&9k9kli - gkgklxSgk+igkli] 



d_ 

~dc 



N-l 



e=0 



k=0 
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N-l N-l 
k=l r=l 

where again we have used discrete integration by parts shifting the sequence gi~ — ► g r 
with r = k + 1, and the fact that <5<7o = 5gjv = 0. Defining dk = Sgkg^ 1 , we obtain 
the discrete Euler-Poincare equations for all variations of this form. □ 



Remark 2.1. In the case thatlL is left invariant, the discrete Euler-Poincare equa- 
tions take the form 

-i'ifkk-i^Rf^ + l'(f k+ i k ) Ad fk+lk TR fk+lk = 0, (2.5) 

where fk+ik = 9k+\9k is in the left quotient (G x G)/G, and the operators act on 
variations of the form — g^Sg^. 

We may associate to any C 1 function F on G x G its Hamiltonian vector field Xp 
satisfying Xp Jwl = dF. The symplectic structure naturally defines a Poisson 
structure {•, -}gxg on G x G by the relation 

{F,H} GxG =uj l (Xf,Xh). (2.6) 

Theorem 2.2. If the action of G on G x G i s pr oper, then the algorithm on G 
defined by the discrete Euler-Poincare equations \2.\ ) preserves the induced Poisson 
structure {•, - }q on G qiven by 

{/, h} G O 7T = {/ O 7T , h o tt} GxG (2.7) 

for any C 1 functions f,h: (G x G)/G = G —> R. 



Proof. Theorem 4.1 of |MPS 98| guarantees that the DEL algorithm preserves the 



symplectic structure wionGx G; hence, by (2.6) ; the DEL algorithm preserves the 



Poisson structure on G x G. Since the action of G on G x G is proper, the general 



Poisson reduction theorem | MR 94 states that the projection 7r:GxG^Gisa 
Poisson map. 



By Theorem 2.1, the projection of the DEL algorithm, 

7T° (Sfc-l,Sfc) >-> 7T° (9k,9k+l) 

is equivalent to the DEP algorithm on G, fk-ik ^ fkk+X- Therefore, as the Poisson 
structure on G is induced by 7r and as 7r is Poisson, we have proven the theorem. □ 

As we shall prove in the following theorem, reconstruction of the DEP algorithm 
(|]J) on G reproduces the DEL algorithm on G x G. 

Theorem 2.3. The discrete Euler- Lagrange algorithm governed by L and the dis- 
crete Euler-Poincare algorithm governed by t are related as follows. The canonical 
projection of a solution of DEL gives a solution of DEP, while the reconstruction 
of a solution of the DEP equations results in a solution of the DEL equations. 

Proof. The first assertion follows by construction. For the second assertion, using 
the definition f kk +i = 9k9k+\, the DEL algorithm can be reconstructed from DEP 
algorithm by 

(9k-i,9k) ^ (9k,9k+i) = (fk-ik • 9k-i,fkk+i ■ 9k), (2.8) 
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where fkk+i is the solution of (2.4). Indeed, / 



'fcfc+i 



each increment, one need only compute f k k+i • gk since gk — f^-ik ' 9k— l is already 



gk is precisely gk+i- Thus, at 
-l 

known. 

Similarly one sho ws that in the case of a left G action, the reconstruction of the 
DEP equations ( |2.5| ) is given by 

(9k-i,9k) >-*■ (gk,gk+i) = {gk-i ■ fkk-n9k ■ fk+ lk )- ( 2 - 9 ) 

□ 



Remark 2.2. Let us denote by W the quotient map W : TG — > TG/G = g mapping 
g € T g G to gg^ 1 G g. In the limit as the time step h — ■> 0, the DEL algorithm 
converges to the flow of the EL equations. 

We denote the reconstruction of the flow of the Euler- Lagrange equations from 
the flow of the Euler- Poincare equations by VKep- Similarly, we denote the recon- 
struction of the DEL algorithm from the DEP algorithm provided by Theorem 2.c 
by y^DEP- The following noncommutative diagram shows these relations. 



GxG 



G 



h— i-O 



TG 







DEL 



DEP 



EL 



EP 



where G x G — * TG as h — > in the following sense. Locally, G x G = iTL* (T* G) 
and as h — > 0, _FL — > FL which pulls-back T*G to TG. Thus, the DEP algorithm 
is asymptotic to the flow of the Euler- Poincare equations if properly interpreted by 
means of reconstruction. 



3. The discrete Lie-Poisson algorithm 

In addition to reconstructing the dynamics on G x G, we may use the coadjoint 
action to form a discrete Lie-Poisson algorithm approximating the dynamics on 
q* . Recall that in the Lie-Poisson reduction setting, for m £ T*G, the momentum 
corresponding to the velocity vector g £ T g G, we define 

m c = L*m G g*, m s = R* g m e g* 

to be the body and spatial momentum vectors, respectively, with the relation 



m s = Ad*_i m c . 

For the right invariant system, the first Euler theorem states that (d/dt)m c = 
(see Theorems 4.4 of Arnold and Khesin | AK 98 |), so that the body momentum is 
a constant of the motion. For convenience, we denote the constant m c by [1q and 
m s (t) by n{t) so that 

H(t) = Ad*-i W -fi - (3.1) 

Now, let O C g be a coadjoint orbit. Then O is a symplectic manifold with unique 
Kirillov-Kostant forms w ± as the coadjoint orbit symplectic structures (see, for 



example, Theorem 14.4.1 in [ |MR 94[| ). Lemma 14.4.2 of |MR 94| states that for 
any j e G, Ad*-i : O — » O preserves w ± . On the other hand, there are natural 
Lie-Poisson {•, •} ± structures on g* (coming from Lie-Poisson reduction on T*G) 
which induce (±) symplectic forms on each symplectic leaf in g*. These induced 
symplectic structures coincide with the coadjoint orbit symplectic structures on 
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each coadjoint orbit (see Kostant [ K 66| |); hence, the coadjoint action preserves the 
Lie-Poisson structures. 



Using the evolution equation (3.1) along with the sequence {fkk+i} obtained by 
the DEP algorithm, we find that 

fi k+1 = Ad*-^ = Ad*^^^., ^ = Ad) kk+1 .Ad^-t mo = Ad} kk+i Mfc . 
Thus, we have proven the following 

Proposition 3.1. An algorithm, called the discrete Lie-Poisson (DLP) algo- 
rithm, on q* defined along the sequence {fkk+i} provided by the DEP algorithm on 
G and given by 

fj,k+i = Ad} fcfc+1 -uk (3.2) 
is Lie-Poisson, i.e. it preserves the (+) Lie-Poisson structure on g* . 

Remark 3.1. The corresponding discrete Lie-Poisson equations for the left invari- 
ant system is given by Q 

n fe+1 = Ad*x -rife, (3.3) 

Jk + lk 

where Uk '■= Ad* fc ttq and the reduced variable m c {t) is denoted by H(t) and the 
constant m s by ttq. 



Thus, one can obtain a Lie-Poisson integrator by solving (2.4) for fkk+i and 
then substituting it into ( |3.2[ ) to generate the algorithm. This algorithm manifestly 
preserves the coadjoint orbits and hence the Poisson structure on g*. In Section 
[|, we shall show that this recovers the Moser-Veselov equations for genenaralized 
rigid-body dynamics on SO(n). 

It is instructive to compare our discrete Lie-Poisson algorithm with that obtained 



by Ge and Marsden [GM 88 1 using the Lie-Poisson Hamilton- Jacobi equations. We 
now state their results which were obtained for the left action of a group G on itself. 
Let H be a G-invariant Hamiltonian on T*G and let Hl be the corresponding 
left reduced Hamiltonian on g* . If a generating function S : G x G — » K of 
canonical transformations is invariant, then there exists a unique function Sl such 
that S^g^go) = S(g,g ). 

The left reduced Hamilton- Jacobi equation for the function Sl ■ G — > R is given 

by 

a o 

-^ + H L (-TR* g -dS L (g)) = 0, (3.4) 

and is called the Lie-Poisson Hamilton- Jacobi equation. The Lie-Poisson flow of 
the Hamiltonian Hl is generated by its solution Sl] in particular, the flow 1 1— > Ft 
of Sl taking initial data Ho to H(i) is Poisson for each t in the domain of definition. 
Next, one defines g G G as the solution of 

n - -TL* g ■ D g S L (3.5) 

and then sets 

n = Ad* 1 n . (3.6) 



1 Henceforth, we shall use the notation fj, £ 0* for the right invariant system and II G g* for 
the left. 
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Thus, one obtains a Lie-Poisson integrator by approximately solving (3.4), and then 
using (3.5) and (3.6) to generate the algorithm. 

Note that (3.4) is the analogue of the usual Hamilton- Jacobi equation 

dS „ ( t dS s 



o 



and that (3.5) and ( p.6[ ) are the analogues of the corresponding canonical transfor- 
mations generated by a solution S which in a local chart are given by 

dS dS 

Oq l Oq 1 



It is interesting to compare the Lie-Poisson Hamilton- Jacobi equation (3.4) with 
the discrete Euler-Poincare equation (2.4). Although it may be less comp uta tionally 
intensive to solve (3.4), one must further solve the difficult equation (3__5) before 
advancing the algorithm by (3.6). Therefore, the discrete Lie-Poisson algorithm 
(3.2) provides an explicit and more efficient method for constructing a Lie-Poisson 
dynamics as compared to the implicit equations (3.5) and (3.6). Namely, although 
both (3.2) and (3.6) manifestly advance the discrete trajectory along the coadjoint 
orbit, the DLP equation provides a time evolution map fj,k i— » Hk+i on q* using a 
known solution fkk+i, while (3.6) advances the initial value Ho along the coadjoint 
orbit and requires at each time step the solution g of (3.5) that approximates the 
current "position" g(t). 

4. Discretization using natural charts 

In this section, we discretize TG by G x G and use the group exponential map 
at the identity, exp e : q — > G, to construct an appropriate discrete Lagrangian. 

4.1. The general theory. For finite dimensional Lie groups G, exp e is locally a 
diffcomorphism and thus provides a natural chart. Namely, there exists an open 
neighborhood U of e G G such that exp^ 1 : U — > u = exp~ 1 ([7) is a C°° diffco- 
morphism (this is not in general true for infinite dimensional groups). Hence, the 
manifold structure is provided by right translation, so that a chart at g 6 G is given 

by 

ipg = exp,: 1 oi? g -i. (4.1) 
We now define the discrete Lagrangian, L:GxG 



L (3i>32) = L V 



»(V) 



■ R, by 
^9(32) 



ipg(9l) 



(4.2) 



where h G R+ is the given time step and 31,32 £ U g = R g (U). 

We shall assume that G has a right invariant Riemannian metric (•, •) obtained 
by right translating a positive bilinear form on q over the entire group. For K C G 
a compact set, we define the Riemannian distance function, dist : K x K — > R + by 

dist( 5l , 52 )= f (j(t),j(t))dt, 
Jo 

where 7 : [0, 1] — > G is the geodesic with 7(0) = g\ and 7(1) = g^. It is then clear 
that diam([/) = diam(J7 g ) for all g G G, so in order for ( |4.2| ) to be well defined 
we require that dist(<?i, g%) < diam(J7). In other words, we require that (31,32) 
be close to the diagonal in G x G. Our restriction on dist (31,32) m turn places a 
restriction on the timestep h. 
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n 



with corresponding group element 

g' = exp(?7) G U. 

We denote the algebra element approximating the velocity g~ l g by 

c _ ^ g (g2) -^g(ffi) 

* h 

Using the standard formula for the derivative of the exponential (see, for example, 
Dragt and Finn | DF 76 or Channel and Scovel CS 91 ) given by 

T v exp = T e R gl ■ iex(- ad r) ), i)£g, g' = exp(r}) G U, 

where iex is the function defined by 



iex(w) = 



- (n+l)V 



(4.3) 



n=0 



we may evaluate the push-forward of V'g 1 at 77. We obtain the following expression 
for the discrete Lagrangian 

Lfoi.sa) = i {^g l {rj),T gl R g ■ T e R g , ■ iex(- ad„)(0) . 

Setting g = ip~ (rj) = R g g\ the last formula is expressed as 



L (si,52) = L(q,T e R q • iex(- ad„)(C)) 



(4.4) 



so that locally the Lagrangian is evaluated at the base point q = ipg 1 ^) G U g C G, 
and the Lie algebra (fiber) element iex(— ad,,)(£) is right translated to the tangent 
space at the point q, T q G; as h — > 0, this fiber element converges to the group 
velocity g £T S G. 

The following lemma establishes that the discrete Lagrangian L inherits the G- 
invariance property from the original Lagrangian L, so that the discrete counterpart 
of the Euler-Poincare reduction is well-defined. 

Lemma 4.1. The discrete Lagrangian L : G x G — > K is right (left) invariant 
under the diagonal action of G on G x G, whenever L : TG — ► K is rig/ii ('ie/tj 
invariant. 

Proof. We fix the right action and consider i?|(L) for some <? G G. By construction, 
R g gi, R s gi G R g (U g ), whenever gi,<?2 G l/ fl = R g (U), so that the chart is given by 

V-Sff = eX Pe 1 °- R ( 9 g)- 1 - 

By defintion, both 77 and £ are always elements of a neighborhood of G g, so 
it is clear that they are right invariant. Hence, using the explicit form of the chart 
ip g g t oge ther with the right invariance of the Lagrangian L, we obtain from (4.2) 
and (fir)) that 



h{R g g 1 ,R g g 2 ) 



i>gg{9l9) +^gg{929) 



> (4>gg) 



^gg{92g) -ipgg(9ig) 



L (R g ■ ^(rj), T q R g ■ T g ,R g ■ T e R g , • iex(- ad r) )(C)) 
L (R g ■ q, T q R g ■ T e R q • iex(- ad„)(0) 

Hgi,g2). 
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In the case that the group action is on the left, we use <f> g 
chart, and proceed with the same argument. 



exp e 1 oL g - 



as the 
□ 



Corollary 4.1. Using the discretization defined by the reduced discrete La- 

grangian I determined by the projection map $2.$) , ^(gig^ 1 ) = M.9ii$2), can be 
expressed in terms of the continuous reduced Lagrangian I by 

£(5i^ 1 )='(iex(-ad^)(0), (4.5) 

where i] = (ip g (gi) + ip g (g 2 ))/2, C = {^ g {9i) - ip g (gi))/h, and I can be defined by 
translation to the identity of the arguments of the right invariant Lagrangian L, i.e. 
1(0 = I Mi,, g,TR g g) = L(e,0, where £ = TR g - ig 6 Q. 



The proof of this corollary follows from expression (4.4), and the fact that the 
Lagrangian L is right invariant so that translation by q~ l to e gives (4.5). 

The expressions (4.4) and (4.5) for the discrete Lagrangian in general require 
evaluation of the infinite series for the iex function given by (4.3); however, a 
simplification occurs when g is set to either gk or gu+i- This is due to the fact that 
when g = gk or g ~ gk+i, one may easily verify that ad^rj := [C, rj\ — 0, and hence 
that iex(— ad^)(C) = C- 

For example, with g = gk+i, the discrete Lagrangian is simply 

Hgk,gk+i) =L(q,T e R q (C)) , 

where 



(4.6) 



V = ^log(3fc5 fe +i), 9 



and log : 



2 
exp 



. Consequently, the reduced discrete Lagrangian is given by 

i{fkk+l) = l(\0g(fkk+l)/h), 

where f k k+\ = 9k9k+i- 



(4.7) 



Substituting the discrete Lagrangian (4/7) into the DEP equation (2.4), we obtain 
the following implicit algorithm on the Lie albebra 



I'Ukk+i/h) ■ x(a.d ikk+1 ) = I'fe-ik/h) ■ x(ad Cf! _ lfc ) • exp(ad 5fc _ lfc ), 



(4.8) 



where £fcfc+i = logfkk+i € fl and the function x is defined to be the inverse of the 
function iex defined by (f4.3j), x(ad^) • iex(— ad^) = Id fl . The function x m (4.8) 
arises from taking the derivative of the log function viewed as a map from the Lie 
group to its algebra. It is interesting to compare the above algorithm with the one 
obtained by Channel and Scovel [CS91] using the Hamilton- Jacobi equation. 



4.2. Generalized rigid body dynamics. We apply our DEP algorithm to the 
generalized rigid body problem. In this case, G — SO(n) with Lie algebra g = so(n), 
and the left invariant Lagrangian is given by the kinetic energy 



1 



1 



1 



1 



(9,9) = ^{9,9)3 = ^{gJM) = 2<<r 1 <?,JGr 1 <?)} = u-u in- (4.9) 

-) g denotes the pairing between T g SO(n) and its dual T* SO(n) which we 



Here, ( 

associate to the metric (•, •) on SO(n) by 

(X g ,Y g ) g = (X g J g Y g ) g , X g ,Y g E T g SO(n), 

where J s = (I/*) -1 J (L g -i)* is the left translated inertia tensor, and J : so(n) 
so(n)*. On SO(n), [L g -i ), ■ g = g^g. 
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We discretize T* SO(n) by SO(n) x SO(n) and construct the discrete Lagrangian 
following ([0]) as 

L_Rs(fffc,5fe+i) = Lrb (qk+ik,T e L qk+lk ((^ k+lk )) , 

where q k+lk = gk+i{gkli9k) 1/2 and ( k+lk = ^ogig^gk)- Using the left invari- 
ance of the metric, we may express the discrete rigid body Lagrangian as 

^RB(gk,gk+i) = 2 (Cfc+ifc: Cfc+ifc) = ^{Ck+ik, J(Cfe+ifc)}- (4.10) 

The Lagrangian for the reduced system on (SO(n) x SO(n))/ SO(n) = SO(n) is 
then given by 

^RB{fk+lk) = L_Rs(.g/c,5fe+l) = ^(log/fc+ifc, J(log/ fe+ i fe )), (4.11) 

where fk+ik = S^+ifffc G SO(n) is an element of the reduced space and h is the 
time step. 



The DEP equation (2J5) has the following implicit form 



Cfc+ifc = J 1 (iex(- ad* hCk+ J ■ x(ad* hCkk i ) • Ad^ p( _ /tffe)s _ l) J(Cfcfc-i)) • (4.12) 

5. Moser-Veselov discretization of the generalized rigid body 

An alternative discretization approach may be taken if we first embed our group 
G into a linear space; for finite dimensional matrix groups, the linear ambient space 
is fll(n). Then, summation of the group elements becomes a legitimate operation 
provided we project the result back onto the group G by using Lagrange multipliers. 

In this section, we consider the left invariant generalized rigid body equations 
on SO(n). The corresponding Lagrangian is determined by a symmetric positive 
definite operator J : so(n) — > so(n), defined by J(£) = A£ + £A, where £ € so(n) 
and A is a diagonal matrix satisfying Aj + Aj > for all i ^ j. The left invariant 
metric on SO(n) is obtained by left translating the bilinear form at e given by 

&0 = jTr(C T J(0). 

The operator J, viewed as a mapping J : so(n) — > so(n)*, has the usual interpre- 
tation of the inertia tensor, and the Aj correspond to the sums of certain principal 
moments of inertia. 

The rigid body Lagrangian is the kinetic energy of the system 

L(g,g) = |<<r 1 <?,J(<r 1 <?)> = (5-i) 

where £ = g~ x g G so(n) and (■, •) is the pairing between the Lie group and its dual; 
hence, the Hamiltonian vector field of L is the geodesic spr ay on TG. 

Using the definition of J we rewrite the Lagrangian ( |5.l|) in the following form: 

We now discretize the Lie algebra elements by £ = g~ Y g 

j^fffe+iOfc+i -5fc)> ( 5 - 2 ) 
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where h is the time step. Substituting ( p.2[ ) into the Lagrangian L (and using prop- 
erties of the trace), we obtain the following expression for the discrete Lagrangian 
(modulo K): 

H9k,9k+i) = -Ta Tr (9k^9k+i) ■ 
We remark that exactly the same expression is obtained if we instead discretize £ 
by Tffk{9k+i ~ 9k)- Notice that up to a multiplier of —1/h 2 , this is precisely the 



Lagranigian used by Moser and Veselov [MoV 91 



We scale the above Lagrangian and introduce matrix Lagrange multipliers Xk, 
imposing the constraint &k{gk) = 9k9k~ Id = 0. By decomposing Xk into symmetric 
and skew components, we see that the skew component of Afc does not contribute 
to the action because the constraint is symmetric; thus, we find that Afc = A^ . 
The action sum then takes the form 

S = M.9fc+i) - \ E Tr (A* {g k gl - Id)) (5.3) 

k k 

Notice that the discrete Lagrangian L is left invariant and can be reduced to 
a Lagrangian i : G — ► K using the canonical projection ir : (ffk>5fe+i) l— * fk+ik = 
9k+i9k so that 

£(fk+ik) = Tr(/ fe+lfe A). 
Because the constraint, ensuring that each gk € G, is G-invariant, there exists a 
Lagrange multiplier Xk in the conjugacy class of Afc, i.e., Aft = g T Xkg for all g G G, 
so that Afc = Afc. Hence, computing the discrete variation of Tr(Afe$fc(pfc)) with 
respect to we obtain the operator equation 

-Z'{fkk-i)TRf kk _ 1 +t'(f k+ i k )Adf k+lh TR fk+lk = X k , 

where the operators act on the variations $k = 9kSgk- Using the expression for the 
reduced Lagrangian t, the DEP equation can then be written as 

/fe+ifeA + /fefe-iA = X k . 

Using the fact that A^ = Afc, we obtain the DEP algorithm on SO(n) as 

fk+ik A - A/fc +lfc = A/ fe r fe _ 1 - /fcfc_iA. (5.4) 

This is an implicit scheme to be solved for fk+ik using the current value fkk-i- 



The solution of (5.4) generates the explicit DLP algorithm on so(n)* given by 



rifc +1 = Ad,-i PL = /fc. +1 fc.nfc/J +u . (5.5) 

•1 k-\- lk 

Finally, reconstruction of the DEP algorithm recovers the DEL algorithm on 



G x G which, according to (2.8), is given by 



{9k-i,9k) i-> (gk,9k+i) = {9k, 9k ■ /fc+ifc)- 



Theorem 5.1. The above DEP and DLP algorithms given by (5.4) and (5.5), re- 
spectively, are equivalent to the Moser- Veselov equations 

A/fc+i = whM^Ij , . 

Mfc = uj^A - Acjfc, uj k £ SO(n), v ■ ' 



where (using the notation of [MoV 91 |J uok = g k 9k-i E SO(n) is the discrete angu- 
lar velocity, Mk = g^— l^fcfffe-i = w fe A — Au>k E so(n) is the discrete body angular 
momentum, and m k = mo is the constant discrete spatial angular momentum. 
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Proof. Comparing the definitions of fkk-x = 9k9k-x an d ^>k — 9k9k-i, we see that 
fkk-i = w fc- Similarly, comparing the definitions of Ilfc = Ad* fc 7Tq and 



M k = g^mkgk-i = 9k-i m o9k- 



Ad! 



m . Hence, the first equation in (|5.6|) is 



we conclude that Uk-i = M k and n Q 
precisely the DLP algorithm ( |5.5| ). 

Substituting the second equation of (5.6) into the first results in the following 
expression: 

w^ +1 A - Auj k+ i = Kul - w fe A, 

which is precisely the DEP equation fl5.4| ) when the above identifications are in- 
voked. □ 

The Moser-Veselov algorithm (^1]) has an an obvious geometric mechanical inter- 
pretation. The first equation can be viewed as a discretization of the left Lie-Poisson 
equation 

M k = 9k~i m o9k-i = Ad* fc l mo, 

rewritten in terms of the uik and this corresponds to the DLP algorithm ( |5.5| ) . The 
second equation is a discrete version of the relation between the angular momentum 
and angular velocity, as it is obtained by substitution of ( |5.2| ) into M = J(£) = 
A£ + £A. 



The DEP algorithm (5.4) provides an equivalent alternative to the Moser-Veselov 
scheme ([5.6|), the difference being that the former is an algorithm on G only, while 
the latter is a combined algortihm on G and q* and schematically can be represented 
by the mappings g* i— » G >— > g* >— > G; Mk i— > i— > M/s+i i— > Wfe+i. 

we identified n^. 



In proof of Theorem 5.1 



_i with Mk in order to establish the 
equivalence with the Moser-Veselov algorithm; however, without any such identifi- 
cation, we exactly obtain the algorithm given by equation (4.1) in Lewis and Simo 
|LS 96| which we write in our notation as 



9k+i — gkfk+ik> 

IIfe+1 = fk+lk^-kfk+lki 

AtUk = 2skew(g fe A). 



(5.7) 



Equation (5.7Q corresponds t o o ur reco nstr uction algorithm (2.S), (p.! 2) corre- 



sponds to our DLP algorithm (3.5), and (5.7 3) is our DEP algorithm (5.4). To see 
this, simply note that 



g{ ([LS 96|, Eq. 4.5) g k = Eq. 5.4 (i.e. DEP) 



It is worthwhile to make a few remarks at this point. Although it is claimed 
in |LS 96 that a computation of the first variation of the action J^ fc Ti(gkAg k+1 ) 
leads to the algorithm (5.7), we have shown that only constrained variations of the 
action function (5.3) lead to this algorithm. Furthermore, the algorithm (5/7) is 
obtained by constraining the iterates of the momentum to be equal; this constraint 
is superfluous as the discrete Euler-Lagrange equations necessarily conserve the 
momentum. Finally, if we choose fk+ik = cay(£fc+ifc) where cay: so(n) — > SO(n) is 
the Cayley transform given by cay(£) = (1 + — f° r an y £ e so( re), then 

the rigid-body algorithm for £k+ik is second-order accurate, as proven in [LS ' 



It is not clear, however, whether the second-order accuracy can be maintained in 
the absence of the Cayley transform. 
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6. A COMPARISON OF DEP/DLP ALGORITHMS WITH SPLITTING METHODS 

For the purpose of comparison, we shall now describe the Hamiltonian splitting 
methods for generating Lie-Poisson integrators on g*, the dual of the Lie algebra of 
a group G. The basic idea behind the construction of such an algorithm follows from 
the fact that many Lie-Poisson systems are governed by reduced Hamiltonians h 
which can be written as a sum h 1 + - ■ - + h N , where each h l can be exactly integrated. 
Letting 4>\ denote the flow of the Hamilonian system h % , we see that to first order 
in the time-step At, the flow <j> t generated by h may be expressed as 

<pAt = 4>A t ° ■ ■ - ° <t>At- 
As each of the maps <fi l At is a Poisson map, hence symplectic on each leaf, the 
composition must also preserve the Poisson structure. Consequently, all Casamirs 
are also preserved by this splitting algorithm. Furthermore, one may construct 
this splitting algorithm to any order of accuracy in At. (For example, the leapfrog 



method 4>±At4>_i At IS a second order accurate scheme (see, for example, [McS 96]).) 

Whereas the DEP/DLP algorithms manifestly preserve the Poisson structure 
and all of the corresponding Casamirs as well, they do much more. First, the re- 
duced algorithms may be used in both the Lagrangian and Hamiltonians sides, in 
that computation of the discrete Euler-Poincare trajectory immediately leads to the 
discrete Lie-Poisson trajectory on g*. More importantly, the discrete Lie-Poisson 
or Euler-Poincare dynamics may be reconstructed to obtain symplectic-momentum 
integrators on TG, for example. Conservation of momentum ensures that the re- 
constructed discrete trajectory lies in an n dimensional submanifold of the full 2n 
dimensional space G x G, approximating TG. This n dimensional submanifold is 
the level set of the discrete momentum mapping. For a small enough time step At, 
G x G is locally diffeomorphic to TG through the discrete Legendre transform, and 
hence we ensure that our discrete reconstruced trajectory is conserving the actual 
momentum. 

Now recall that for right invariant systems, we have used the variable m s to de- 
note the solution of the Lie-Poisson equation, from which we obtain that m c (t) = 
Ad*( ( ) m s (t) is conserved. Using our DEP algorithm, we may compute the dis- 
crete trajectory {(m s ) kk+1 } , reconstruct to find gk, and find that (m c ) tj!+1 = 
Ad* fc (m s ) kk+1 is conserved. On the other hand, the splitting method does not pro- 
vide an algorithm for reconstructing the motion on T*G in such a way as to ensure 
conservation of momentum; thus, there is no obvious way to define the discrete 
analogue of m c , let alone check that it is conserved. 

Nevertheless, there are some computational advantages to using the splitting 
method; the fact that the splitting method leads to an explicit scheme is perhaps the 
most important of these advantages. An efficient excplicit algortihm for the SU(n) 
model of two dimensional hydrodynamics on a torus is constructed in [ Mc 93| . The 



author presents a Poisson integrator of complexity 0(N 3 logiV) which preserves 
N — 1 Casimirs. 

7. Addendum: relation to other works 

It is very interesting to compare the above constructions and algorithms to the 
recent results of Bobenko and Suris |BS 98 . In this paper they consider the theory 



of discrete time Lagrangian mechanics on Lie groups and, more specifically, address 
the issue of discrete Lagrangian reduction using left or right trivializations of the 
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(co)tangent bundles of Lie groups. They adopt a somehow broader point of view 
when the symmetry group of a system defined on a Lie group G is a subgroup 
of G. Hence, it includes the Lie-Poisson case as a special case. Below we shall 



demonstrate that the reduced discrete equations obtained in [ |BS 98| agree with 
our DEP/DLP algorithms when the symmetry group is taken to be the full group 
G. Here we summarize their results choosing for consistency and simplicity the 



case of right trivialization and refer the reader to BS 98] for details of proofs and 
notations. 

Let the discrete Lagrangian *L(g k , g k +i) '■ G x G — > K define a a discrete system 
with the corresponding DEL equations. Consider the map 

(g k ,w k ) EG x G i-> (g k ,g k+ i) E G x G, (7.1) 

where 

g k+ i = w k g k <^> w k = g k +\g k l . 
Consider also the right trivialization of the cotangetn bundle T*G: 

(g k ,m k ) EG x * i > (g k ,U k )ET*G, 

where 

n fc = R* im k 44> m k = R* k ir k . 



Denote the pull-back of the Lagrange function under (7.1) by 

lS r) {g k ,w k ) = h(g k ,g k+1 ). 



Proposition 3.5 of BS 98| gives the DEL equations in these coordinates: 

f A ^w k m k+i = m k + d g lS r Xg k ,w k ), < 72 , 
\ g k+ i = w k g k , 

where 

m k = ^L (r) (j fc _i,ii)n) G q*. 
Assume that for some ( £ g, IL.W is invariant under the action of a subgroup 
C?K1 = {h E G\kd h C, = (} C G on G x G induced by right translations on G: 

li r \gh,w) =L^(g,w), hEG [c] . 
Define the reduced Lagrange function A( r ) : x G i— > K. as 

A (r) (a,u>) =L (r) (. 9 ,u>), a = Ad 9 Ceg C 
here 2c is the adjoint orbit of £. 



Then, proposistion 3.7 of [BS 98| states that under the reduction by G^l, the 



reduced Euler-Lagrange equations become 



A C fc m fe+i = m k - ad^V Q A( r )(a fe , w k ), ^ ^ 

a k +i = A.d Wk a k , 



where 

m k = dwA^fak-uWk-i) G 0* 



In (7.3) the follwoing notations are used (see [BS 98|). For a function / : G 



its left and right Lie derivatives, df(g) : G i— > 0* and d'f(g) : G i— » 0*, are defined 
via 

(df(g),v) = jj{^9)\^ V r, E 0, 
(d'f(g),r,) = ±f(ge^)\ e=0 , y V E 5 . 
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Then, the gradiaent V/ : G h-> T*G is related to the above derivatives via 

Vf(g) = R;- 1 df(g)=L* g ^d'f(g). 

Notice that in the Lie-Poisson case, when the symetry group is G itself, the 
r edu ced space is simply the group G represented by w k = g k +ig k l , and equations 



(7.5) become 

M* Wk m k+ i=m k (7.4) 

with 

m k = dA^K-i) G 0*. (7.5) 

Comparing the above notations with the results in our paper, we immediately see 
that w k correspond to the other choice for the quotient map ( |2.3[ ) tt : (g k ,g k +i) l— * 
A-fe+i = 9k9k+i> i- e - fkk+i = Wfc 1 . Similarly, the reduced Lagrangian A^iwk) 
corresponds to l(fkk+i) i n °ur notations. Finally, using the definitions of the Lie 



derivatives above we obtain for the angular momentum (7.5) 



m k = dA«(t«fc_i) = iC^VAM^-i) = R* 1 £'(f k - lk ), 

J k 



where we have substituted our notations. Hence, (7.4) can be written as 



Ad* x IT x £'(fkk+i) = fl* r i ^(A-ifc). 

Jfcfc+l Jfcfc+l Jk-lk 



The last expression is precisely the DEP algorithm fl2.4|) after rewriting it with 
the adjoints of the above operators acting on the variation i9 k = 5g k g k (see section 



2). It is interestig to note that the second equation in (7.2) corresponds to our 



reconstruction equation (2.8). Similar correspondence can be estableshed for the 



case of left trivialization considered in [BS 98 
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